import os
import numpy as np
import bigfish
import bigfish.stack as stack
import bigfish.detection as detection
import bigfish.multistack as multistack
import bigfish.plot as plot
print("Big-FISH version: {0}".format(bigfish.__version__))
Big-FISH version: 0.6.2
# hard-code the paths of our input and output directories
path_input = "../data/input/DMSO"
#path_input = "../data/input/JQ1"
#path_input = "../data/input/TSA"
path_output = "../data/output/DMSO"
# path_output = "../data/output/JQ1"
# path_output = "../data/output/TSA"
In this notebook, we show examples of mRNAs detection. We use 3D images. Three main steps are developed:
#Change filename
path = os.path.join(path_input, "DMSO_004_smfish_fov_1.tif")
#path = os.path.join(path_input, "JQ1_317_smfish_fov_1.tif")
#path = os.path.join(path_input, "TSA_615_smfish_fov_1.tif")
rna = stack.read_image(path)
print("smfish channel")
print("\r shape: {0}".format(rna.shape))
print("\r dtype: {0}".format(rna.dtype))
smfish channel shape: (15, 512, 512) dtype: uint16
rna_mip = stack.maximum_projection(rna)
print("smfish channel (2D maximum projection)")
print("\r shape: {0}".format(rna_mip.shape))
print("\r dtype: {0}".format(rna_mip.dtype))
smfish channel (2D maximum projection) shape: (512, 512) dtype: uint16
We assume spot is a local maximum in the smFISH channel. Three steps are required to detect them:
bigfish.stack.log_filter).bigfish.detection.local_maximum_detection).bigfish.detection.spots_thresholding). To be robust, the thresholding should be applied on the filtered image. Thus, the threshold is set relatively to the filtered image values.bigfish.detection.automated_threshold_setting (applied on a filtered image).All these steps are summarized in bigfish.detection.detect_spots that return the 2D or 3D coordinates of the detected spots.
spots, threshold = detection.detect_spots(
images=rna,
return_threshold=True,
voxel_size=(300, 103, 103), # in nanometer (one value per dimension zyx)
spot_radius=(350, 150, 150)) # in nanometer (one value per dimension zyx)
print("detected spots")
print("\r shape: {0}".format(spots.shape))
print("\r dtype: {0}".format(spots.dtype))
print("\r threshold: {0}".format(threshold))
detected spots shape: (768, 3) dtype: int64 threshold: 77.0
Given the voxel size and the expected spot radius (in nanometer), the function bigfish.detection.detect_spots automatically estimates a kernel size for the LoG filtering and a minimal distance between two spots we want to be able to detect separately. It is still possible to set these parameters explicitly in order to fine-tune the detection. Internally, we approximate them as the spot radius in pixel with the function bigfish.detection.get_object_radius_pixel.
spot_radius_px = detection.get_object_radius_pixel(
voxel_size_nm=(300, 103, 103),
object_radius_nm=(350, 150, 150),
ndim=3)
print("spot radius (z axis): {:0.3f} pixels".format(spot_radius_px[0]))
print("spot radius (yx plan): {:0.3f} pixels".format(spot_radius_px[-1]))
spot radius (z axis): 1.167 pixels spot radius (yx plan): 1.456 pixels
spots, threshold = detection.detect_spots(
images=rna,
return_threshold=True,
log_kernel_size=(1.167, 1.456, 1.456),
minimum_distance=(1.167, 1.456, 1.456))
print("detected spots")
print("\r shape: {0}".format(spots.shape))
print("\r dtype: {0}".format(spots.dtype))
print("\r threshold: {0}".format(threshold))
detected spots shape: (770, 3) dtype: int64 threshold: 77.0
Note: What we call spot radius in this notebook can be understood as its Point Spread Function (PSF). For simplicity sake, this PSF is modelled as a 2D or 3D gaussian.
plot.plot_detection(rna_mip, spots, contrast=True)
The automated spot detection method tries to find the optimal threshold to discriminate actual spots from noisy blobs. If we plot the number of the spots detected as a function of threshold level we observe an elbow curve. The selected threhold is the one located in the breaking point of the curve. This curve can be plotted with bigfish.plot.plot_elbow.
plot.plot_elbow(
images=rna,
voxel_size=(300, 103, 103),
spot_radius=(350, 150, 150))
The detection of local maximum is not able to detect individual spots clustered in a dense and bright region. We try to decompose these regions by simulating as many spots as possible until we match the original region intensity. Our current steps are:
bigfish.stack.remove_background_gaussian).bigfish.detection.build_reference_spot).bigfish.detection.modelize_spot).bigfish.detection.get_dense_region).bigfish.detection.simulate_gaussian_mixture).All these steps are summarized in bigfish.detection.decompose_dense that return the 2D or 3D coordinates of the detected spots outside and inside a decomposed region, additional information about the regions themself and an image of the reference spot estimated.
spots_post_decomposition, dense_regions, reference_spot = detection.decompose_dense(
image=rna,
spots=spots,
voxel_size=(300, 103, 103),
spot_radius=(350, 150, 150),
alpha=0.7, # alpha impacts the number of spots per candidate region
beta=1, # beta impacts the number of candidate regions to decompose
gamma=5) # gamma the filtering step to denoise the image
print("detected spots before decomposition")
print("\r shape: {0}".format(spots.shape))
print("\r dtype: {0}".format(spots.dtype), "\n")
print("detected spots after decomposition")
print("\r shape: {0}".format(spots_post_decomposition.shape))
print("\r dtype: {0}".format(spots_post_decomposition.dtype))
detected spots before decomposition shape: (770, 3) dtype: int64 detected spots after decomposition shape: (1047, 3) dtype: int64
plot.plot_detection(rna_mip, spots_post_decomposition, contrast=True)
Two spots are considered connected if they localized within a specific radius (in nanometer). Above a minimum number of connected spots, a cluster can be defined. This detection can be computed with bigfish.detection.detect_clusters.
spots_post_clustering, clusters = detection.detect_clusters(
spots=spots_post_decomposition,
voxel_size=(300, 103, 103),
radius=350,
nb_min_spots=2)
print("detected spots after clustering")
print("\r shape: {0}".format(spots_post_clustering.shape))
print("\r dtype: {0}".format(spots_post_clustering.dtype), "\n")
print("detected clusters")
print("\r shape: {0}".format(clusters.shape))
print("\r dtype: {0}".format(clusters.dtype))
detected spots after clustering shape: (1047, 4) dtype: int64 detected clusters shape: (113, 5) dtype: int64
# plot
plot.plot_detection(rna_mip,
spots=[spots_post_decomposition, clusters[:, :3]],
shape=["circle", "polygon"],
radius=[3, 6],
color=["red", "blue"],
linewidth=[1, 2],
fill=[False, True],
contrast=True)
Spots and foci coordinates can be saved in npy files (numpy dedicated format) or csv files using functions bigfish.stack.save_array and bigfish.stack.save_data_to_csv respectively.
# save in csv files
path_spots = os.path.join(path_output, "DMSO_004_spots.csv")
#path_spots = os.path.join(path_output, "JQ1_317_spots.csv")
#path_spots = os.path.join(path_output, "TSA_615_spots.csv")
stack.save_data_to_csv(spots_post_clustering, path_spots)
path_clusters = os.path.join(path_output, "DMSO_004_clusters.csv")
#path_clusters = os.path.join(path_output, "JQ1_317_clusters.csv")
#path_clusters = os.path.join(path_output, "TSA_615_clusters.csv")
stack.save_data_to_csv(clusters, path_clusters)